# rmarkdown::render("./Scripts/03_explore_urban.R", output_dir = "./Results/")
library(here)
library(tidyverse)
library(sf)
## Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
library(rnaturalearth)
library(fasterize)
##
## Attaching package: 'fasterize'
## The following object is masked from 'package:graphics':
##
## plot
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:tidyr':
##
## extract
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(mapview)
zoom <- data.frame(lon = c(-80, -75, -80, -75),
lat = c(35, 35, 40, 40)) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(3347) %>%
st_bbox()
na <- ne_countries(continent = "North America", returnclass = "sf") %>%
st_transform(3347) %>%
st_crop(zoom)
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries
urban <- ne_load(scale = 50, "urban_areas", destdir = here("./Data/Raw/ne_urban/"), returnclass = "sf") %>%
st_transform(3347) %>%
st_crop(zoom) %>%
st_cast("MULTIPOLYGON")
## OGR data source with driver: ESRI Shapefile
## Source: "/home/steffi/Projects/Urban Birds/Data/Raw/ne_urban", layer: "ne_50m_urban_areas"
## with 2143 features
## It has 4 fields
## Integer64 fields read as strings: scalerank
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries
ggplot(na) +
geom_sf() +
geom_sf(data = urban, colour = "orange", fill = "grey")
r <- fasterize::raster(urban, res = 500)
urban_raster <- fasterize(urban, r, fun = "max")
urban_plot <- as.data.frame(urban_raster, xy = TRUE) %>%
filter(!is.na(layer))
ggplot(na) +
geom_sf() +
geom_tile(data = urban_plot, aes(x = x, y = y), colour = "orange")
zoom <- data.frame(lon = c(-79.75, -79.5, -79.75, -79.5),
lat = c(43.75, 43.75, 44, 44)) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(3347) %>%
st_bbox()
Not too bad a rasterization job!
ggplot(na) +
geom_sf() +
geom_sf(data = urban, colour = "red", fill = "grey") +
geom_tile(data = urban_plot, aes(x = x, y = y), colour = "orange", alpha = 0.8) +
coord_sf(xlim = zoom[c(1,3)], ylim = zoom[c(2,4)])
Is this 500m resolution? Hmm, approximately (may depend on where around the world)
mapview(urban_raster, maxpixels = 1393639)
Perhaps worth requesting original?